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ABSTRACT 


Nearshore depths for Waimanalo Beach, HI, are extracted from optical imagery, taken by 
the WorldView-2 satellite on 31 March 2011, by means of automated wave kinematics 
bathymetry (WKB). Two sets of three sequential images taken at intervals of about 10 
seconds are used for the analyses herein. Water depths are calculated using a computer 
program that registers the images, estimates the currents, and then uses the linear 
dispersion relationship for surface gravity waves to estimate depth. Depths are generated 
from close to shore out to about 20 meters depth. Comparisons with SHOALS LIDAR 
bathymetry values show WKB depths are accurate to about half a meter, with R 2 values 
of 90%, and are frequently in the range of 10-20 percent relative error for depths ranging 
from 2-16 meters. 
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I. INTRODUCTION 


A. PURPOSE OF RESEARCH 

Characterizing the environmental parameters of the battlespace ahead of a 
military operation can enable the planners and decision makers to understand the 
potential impacts the environment might have. This can help them conduct safer and 
more efficient operations. One parameter that is extremely important for certain 
operations conducted in the nearshore is bathymetry. Amphibious landings, mine 
warfare operations, reconnaissance missions, and other special operations missions 
performed in the nearshore region require accurate and up-to-date bathymetric 
information. Without it, the operation could be hampered by difficulties, such as vessel 
groundings, causing equipment damage, personnel injuries or death, delays, and perhaps 
mission failure. 

Several methods and techniques are currently employed to detennine nearshore 
bathymetry. Hydrographic surveys, especially those that require International 
Hydrographic Organization (IHO) standards, such as for navigation charts, are typically 
performed using acoustic systems, such as single and multi-beam sonar, mounted on a 
vessel. Relying on relatively high frequency sound to illuminate the bottom, they are 
very accurate but time consuming owing to their small swath, sometimes taking many 
days to map a region of interest. 

The Compact Hydrographic Airborne Rapid Total Survey (CHARTS) system 
(Figure 1), an airplane with a Light Detection and Ranging (LIDAR) instrument mounted 
on it, is also very accurate and can perform large area surveys in a much shorter time. 
Both of these methods, however, require uncontested access to the region of interest. 
Often times, military operations are conducted in denied areas, where survey vessels and 
aircraft are in danger of being fired upon. In addition, due to the high cost and demand of 
these scarce resources, it is impractical to resurvey the same area regularly using these 
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methods. This leads to significant time intervals between surveys allowing natural 
nearshore processes, such as tidal currents and storms, to alter the bathymetry, potentially 
rendering older surveys inaccurate. 
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Figure 1. Compact Hydrographic Airborne Rapid Total Survey (From Joint 
Airborne LIDAR Bathymetry Technical Center of Expertise, 2011). 
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If IHO standards are not required, other methods exist that are better protected 
from adversaries through greater stand-off distances, are much lower in cost, and have 
frequent revisit capability, especially those that use satellite-based sensors. Nearshore 
bathymetry can be extracted from satellite multispectral imagery by applying the linear 
dispersion theory of surface gravity waves, which is the focus of this research. Previous 
research into this method at the Naval Postgraduate School focused on determining the 
viability and accuracy of this method, but used a manual, time-intensive process (Myrick, 
2011; McCarthy, 2010). The purpose herein is to take the next step by investigating an 
algorithm that automates this process. A few of the other methods that have been or are 
still being explored are briefly discussed in Chapter II. 

B. SPECIFIC OBJECTIVES 

The objective is to investigate an automated method for determining nearshore 
bathymetry using remotely sensed images of the Waimanalo Beach area of Hawaii. A 
desire to reduce the time and effort required to extract the bathymetry infonnation, 
therefore greatly increasing the efficiency of the process as compared to the previous 
work is what provided motivation for this effort. 

The WorldView-2 satellite took the multispectral and panchromatic images in 
rapid succession in March 2011. The imagery was orthorectified, converted to Universal 
Transverse Mercator (UTM) coordinates, and processed using a Wave Kinematics 
Bathymetry (WKB) computer algorithm to: crop the images to the desired analysis area; 
register the images to one another (~ 1 m accuracy); apply filters and distinguish between 
land, clouds, and water and; finally, extract current and depth fields by fitting a solution 
to the data (Abileah, 2006). The resulting estimated depth fields were compared to 
ground truth data collected by the Joint Airborne LIDAR Bathymetry Technical Center of 
Expertise (JALBTCX) during survey operations using the CHARTS system with a 
Scanning Hydrographic Operational Airborne LIDAR Survey (SHOALS)-3000 
bathymetric LIDAR (figure 1). This study reaffirms that applying the linear dispersion 
relation for surface gravity waves to wave data collected from multispectral satellite 
imagery is a viable technique for determining bathymetry in denied or restricted areas. 
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II. BACKGROUND 


A. INTRODUCTION 

For many years, nearshore bathymetry has been detennined without the use of 
direct measurement, such as a sounding line, even as far back as World War I (Myrick, 
2011). Described in more detail in Myrick’s paper, those early methods consisted of the 
waterline, transparency, wave celerity, and wave period methods. The waterline method 
relied on aerial photographs of the beach at times of low tides, when the bottom was 
partially exposed. The transparency method exploited the concept that deeper water 
absorbs more light and, thus appears darker in aerial photographs. The wave period 
method used the fact that the period is constant as a wave propagates so a relationship can 
be established between the period, the wavelength, and the water depth at multiple 
locations. These methods all had limitations that made them marginally useful, but were 
sometimes arguably better than having no bathymetry information at all (Myrick, 2011). 

The wave celerity method, which invokes the linear dispersion relation discussed 
shortly, was also very limited before more accurate timekeeping and better resolution 
images became available. Modern remote sensing technology has removed these 
previous handicaps and the wave celerity method is now a viable one for detennining 
water depth in nearshore regions (Myrick, 2011). The wave celerity method is the 
method used as the foundation for this study. 

B. THEORY 

Surface gravity waves propagating in the ocean obey the linear dispersion relation 
between wave celerity, or phase speed, c; wave period (7); wavelength (/.); and water 
depth (d). The dispersion relation for surface gravity waves is: 

co 1 = gk tanh(M), (1) 

where co is the angular frequency (co = 2k/I), g is the acceleration due to gravity, and k is 
the wavenumber (k = 2x/X) (Herbers, 2003). For large kd, which is the case in deep 
water, tanh (kd) ~ 1 and the dispersion relation reduces to: 
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CO 2 =gk. 


(2) 


In shallow water, where kd« 1, tanh (kd) ~ kd, so the dispersion relation reduces to: 

co-yfgdk . (3) 

The corresponding limits of the phase speed c = co/k are: 


CO 

for deep water 

(4) 

cs 

II 

toTl 

^4 1 

for shallow water. 

(5) 


The conclusion from Equations (4) and (5) is that deep water waves are 
dispersive, whereas shallow water waves are nondispersive. That is, the phase speed of 
deep water waves is a function only of co, whereas the phase speed of shallow water 
waves is independent of co and is instead, only a function of depth (Herbers, 2003). 
Another way to think of this is, as waves travel from deep to shallow water, their 
frequency and period remain constant, forcing phase speed to decrease and wavenumber 
to increase with decreasing depth. Since these changes are proportional to each other, 
this is exactly the phenomenon that is exploited to extract the water depth. 

In addition to waves slowing and getting shorter as they shoal, another physical 
effect is an increase in their amplitude. This effect becomes important in the surf zone, 
where the waves break and wave height increases wave speed. The nonlinear processes 
that occur in the surf zone and add speed to the waves are not accounted for in the linear 
dispersion relation. This introduces error, usually in the form of overestimated depths 
because of the additional speed (Myrick, 2011). Equation (1) is actually a very good 
approximation for depths greater than 2 m, and is still valid at 1 m depths with moderate 
wave heights (Abileah, 2006). Therefore, caution is necessary when using any method 
that relies on the dispersion relation in the very shallow depths of the surf zone. 

It can be shown that the wave induced velocity components and the water 
pressure decay exponentially with depth as kd —* oo (deep water). In fact, by only half a 
wavelength below the surface, these parameters are reduced to about 4% of their surface 
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values (Herbers, 2003). For this reason, a good thumb rule is: regions where the depth is 
greater than half the wavelength are deep water and the dispersion relation does not help 
in determining water depth, since phase speed does not depend on depth in this regime. 
For the purposes of this study, regions where the depth is less than half the wavelength 
will be considered “nearshore” and is the regime where the dispersion relation is useful 
for this technique. It should be noted that this definition includes both the shallow water 
regime, discussed previously, and the intennediate water depth regime, where phase 
speed depends on both frequency and depth. 

C. SOME MODERN METHODS FOR DETERMINING NEARSHORE 

BATHYMETRY 

Several methods using different techniques have been developed that use modern 
remote sensing technology to detennine nearshore bathymetry. The methods discussed 
are only a few of the many methods and techniques currently being investigated by a 
multitude of researchers. 

1. Airborne Passive Optical System (Using Linear Dispersion Theory) 

The basis of modem methods using the dispersion relation is to collect ocean 
images, extract the space-time characteristics of the waves in the images, then transform 
this data into spectra that can be used to retrieve depth and currents. 

The first method employs a passive optical system mounted on an aircraft (Figure 
2). This turret-based system, called the Airborne Remote Optical Spotlight System 
(AROSS), maps a time series of images to a common geodetic surface by carefully 
measuring the imaging geometry (Dugan et al., 2001a). It was designed using 
commercial off-the-shelf technology with the intention of ultimately mounting it on 
unmanned aerial vehicles. 

Precisely registering the images to each other is critical to obtaining accurate 
results. Other important considerations are: adequate spatial resolution to resolve the 
smallest waves (the resulting viewing geometry typically yields a sufficient 2 m), large 
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enough field of view to capture several of the longest waves (2 km x 2 km is typically 
used), and dwell times long enough to observe several of the longer wave periods (30 s or 
greater is used) (Dugan et al., 2001a). 



Figure 2. Close-up photo of camera and turret as mounted on the nose cone 
of the Pelican aircraft (From Dugan et ah, 2001a). 

One data set using AROSS was collected during the Shoaling Waves Experiment 
in the Outer Banks of North Carolina in 1999. A three-dimensional (3-D) Fourier 
transfonn was used to turn a data cube consisting of 2 minutes’ worth of data from an 
area approximately 500 x 500 m into a frequency-wavenumber power spectrum. Figure 3 
shows a two-dimensional (2-D) slice through the resulting power spectrum, oriented in 
the direction of the primary swell waves. The theoretical dispersion surface that includes 
no current and is for infinite depth is represented by the dashed curve at the intersection 
of the dispersion surface with the plane of the slice. The wave energy is concentrated in a 
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narrow ridge that is slightly skewed off of the theoretical curve. The depth and current 
are detennined by calculating a surface that fits the measured wave spectrum. The 
resulting dispersion surface is denoted by the solid curve and, in this case, a depth of 
8.2 m and current speed of 0.14 m/s is the solution (Dugan, Piotrowski, & Williams, 
2001b). 
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Figure 3. Frequency-wavenumber slice, from Outer Banks study, through 

the power spectrum, oriented in the direction of the primary swell waves. 
The theoretical dispersion surface that includes no current and is for infinite 
depth is represented by the dashed curve. The depth and current are 
determined by calculating a dispersion surface that fits the measured wave 
spectrum, denoted by the solid curve (From Dugan et ah, 2001b). 


Another AROSS study was conducted in the Monterey Bay in 1999. The same 
size tiles were used, but in this case, only 1 minute of data was Fourier transformed. A 
slice through the 3-D spectrum in the direction of the wind waves is shown in Figure 4. 
Again, the solid curve is the intersection of the deep water dispersion surface (i.e., no 
current and infinite depth), but this time the wave energy is closely distributed along this 
curve. Since the data closely match the theoretical curve, it can be surmised that no 
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significant surface currents were present during the collection. There is, however, a 
slight shift of the low frequency energy off the curve. This is due to the longer waves 
(50-100 m wavelengths) feeling the bottom because the depth at the imaged location is 
not infinite (i.e., is less than half a wavelength deep). The actual depth, as detennined 
from nautical charts, is about 18 m (Dugan et ah, 2001a). 



Wavenumber (cycles/m) 

Figure 4. Frequency-wavenumber slice, from Monterey Bay study, through 
the power spectrum, in the direction of the wind. The solid curve is the 
intersection of the deep water dispersion surface (i.e., no current and infinite 
depth). The wave energy is closely distributed along this curve, indicating 
that no significant surface currents were present during the collection (From 

Dugan et ah, 2001a). 


One of the conclusions from these and similar studies is that a dwell (total image 
time series length) of 1 minute or more is required for high accuracy. The root mean 
squared (RMS) error approaches 5% for 2-minute dwell times, but increases rapidly with 
shorter dwell times (Figure 5). This graph is based on using 3-D Fourier transforms in 
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the inversion algorithm. The algorithm in this study uses a 2-D Fourier transform that 
results in a much flatter curve where errors do not sharply increase at shorter dwell times. 
In the aforementioned two studies, the depths agreed within 15% and the currents within 
10% of the in situ values (Dugan, Piotrowski, & Williams, 2002). 
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Figure 5. RMS error as a function of dwell (From Dugan et ah, 2002). 

2. Synthetic Aperture Radar (Using Wave Refraction) 

Another remote sensing application for determining nearshore bathymetry uses 
Synthetic Aperture Radar (SAR), which can be done from both airborne and space-based 
platforms. The basis for this method is to examine wave refraction as determined by 
SAR imagery using the two-scale Bragg scattering model. The Bragg scattering model 
assumes that the SAR image brightness of a patch of ocean is proportional to the 
amplitude of the Bragg waves. Bragg waves are small-scale ocean surface waves that 
have wavelengths equal to that of the projection of the transmitted radar electromagnetic 
wavelength onto the surface, and are propagating directly toward or away from the sensor 
(Wackerman & Clemente-Colon, 2004). An expression for this geometry is: 

A b =-^—, 

2 sin# 
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where X B is the Bragg wavelength, X e is the SAR electromagnetic wavelength, and 9 is the 
local incident angle of the ocean surface. This geometry produces a scattering pattern of 
the SAR radiation that leads to constructive interference and the corresponding brightness 
in the image (Wackerman & Clemente-Colon, 2004). 

The two-scale aspect of the model is illustrated in Figure 6, which shows a slice 
through a simplified ocean surface consisting of one large-scale and one small-scale 
wave. The Bragg waves are the smaller scale waves shown embedded in and riding on 
the large-scale surface waves. Flat plates, a few Bragg wavelengths in size, are then used 
to model the radar’s interaction with the surface. The plates tilt and move based on the 
local slope and motion of the large-scale wave surface, changing the radar image 
signature as they do (Wackerman & Clemente-Colon, 2004). 



Figure 6. Two-scale model illustrating Bragg waves (the smaller scale 
waves) embedded in and riding on the large-scale surface waves. Flat 
plates, a few Bragg wavelengths in size, model the radar’s interaction with 
the surface. The plates tilt and move based on the local slope and motion of 
the large-scale wave surface (From Wackerman & Clemente-Colon, 2004). 


12 
















Bathymetry is estimated using SAR by detennining the amount of wave refraction 
that occurs in the images and using propagation models that predict wave refraction as a 
function of depth. Waves refract as they move from deep water toward the shore at an 
angle to the bathymetric contours, which forces them to turn and align with the shore. 
This turning is due to a differential interaction with the ocean bottom (i.e., one end of the 
wave front reaches shallower water before the other end and, therefore, slows down 
sooner). The faster the depth decreases, the more quickly the wave will turn and become 
parallel to the coast. Since refraction will not occur if the waves start out parallel to the 
shore in deep water, this method is ineffective in those instances (Wackerman & 
Clemente-Colon, 2004). 


3. Aerial Photography (Using Radiometric Techniques) 

A third modem method relies on radiometric techniques to extract water depth in 
shallow areas from multispectral optical imagery. These techniques are based on 
applying algorithms to interpret the radiance received by the imager. This method 
requires fairly clear water to work, so it is not useful in the many turbid nearshore 
locations around the globe. 

A simple water reflectance model that accounts for most of the received signal is 
represented by: 


L i = L SI + kr Bj e 


~Kjfz 


( 7 ) 


where L si is the radiance over deep water; k, is a constant that accounts for the solar 
irradiance, the transmittance of the atmosphere and the water surface, and refraction at 
the surface; r«, is the bottom reflectance; if is the effective attenuation coefficient for 
water; / is a geometric factor to account for the path length through the water; and z is the 
water depth (Lyzenga, 1978). Using this equation to solve for z is the obvious way to 
determine the depth. 

A ratio algorithm was developed to help remove the effect of changes in bottom 
reflectance on the depth calculation, so the model can be applied to scenes in which the 
bottom composition changes. The algorithm can be represented by: 
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where R is the ratio of the bottom-reflected signals in two bands: 

R = (L l -L sl )/(L 2 -L s2 ), 


( 8 ) 

(9) 


and Rb is the ratio of the bottom reflectances in the same two bands. The assumption is 
made that the two bands chosen have the same bottom reflectance ratio for all the bottom 
types in the scene, which is a decent assumption for multispectral systems. This 
algorithm is used with moderate success over relatively clear water on data collected 
from both satellites and aircraft (Lyzenga, 1978). 

To increase operational flexibility and improve performance, a more general, 
albeit more complex, algorithm was developed that modified the preceding model to 
include the effects of scattering in the water and internal reflection at the water surface. 
It depends on proper classification of the bottom types in the scene, which is done prior 
to the depth computation, and it is at least as good as the ratio model at doing this for the 
case of two bands (Lyzenga, 1978). Arguably, using more bands would improve the 
performance of the bottom classification portion of the algorithm. 



Figure 7. Depth error comparison of the ratio and modified methods 

(From Lyzenga, 1978). 
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Once the bottom is classified, the depth, z, is calculated using: 

Y n = B m - Cz , (10) 

where Yn is the set of depth dependent variables, B m is a function of the bottom 
composition (and its value is detennined via the aforementioned classification), and C 
characterizes the water attenuation. This procedure introduces error both through noise 
associated with Y N and if there happens to be any misclassification of the bottom 
(Lyzenga, 1978). Compared to the ratio method, which is subject to much higher noise 
errors, its perfonnance is much improved (Figure 7). 
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III. EXPERIMENT DESIGN 


A. PROBLEM DEFINITION 

Using satellite optical imagery to view wave motion in order to extract the 
bathymetry of nearshore areas using the linear dispersion relation for surface gravity 
waves has its advantages in certain situations. When IHO standards are not required, it is 
much more cost effective and efficient to use this method rather than LIDAR-equipped 
aircraft or sonar-equipped vessels. There is no need to dispatch expensive, hard-to-get 
assets to remote locations, since the existing low earth orbiting satellites will be overhead 
most locations every few days. It is especially good for denied or hostile areas because 
the satellites do not attract attention and remain safely out of the range of most weapons, 
unlike aircraft and ships. 


Table 1. Comparison of bathymetry methods using remote sensing 

(From Abileah, 2006). 



Depth 

accuracy 

Horizontal 

resolution 

Airborne/Satellite 

Turnaround 
(mouths ) 1 

Cost / 
sq kin 1 

Turbid 

waters 

Day/ 

night 

Airborne 

IHO-l 

IHO-l 

Airborne 

6-12 

>$100 

X 

✓ 

Bathymetric Lidar 





0 



Photobathymetry 

Variable 

10-30 m 

Both 

1-2 

$50 

X 

X 

Currents interaction 

Variable 

25 nr 

Both 

1 

$10 

✓ 

V 

Wave refraction 

10-20% 

1 km 

Both 

1 

$10 

✓ 

V 

WKB 

5-10% 

50-200 m 

Both 

1-2 

$50 

✓ 



Notes: Aough-ordei-of-magntude assuming satellites, except ABL. which is on aircraft. ~Ref: Argoss website (mTO~.argoss. 11 l 1. 

3 WKB with the satellite imagery described in this paper is daytime only. However nighttime is possible with long wave IR (Farber. 
1995). nautical X-band radars (see text), and theoretically possibly with space SAR. 


Other satellite remote sensing techniques offer these same advantages. However, 
they each have limitations that make each method suitable only in certain situations. 
Table 1 compares some of the properties of each method. The optical radiance method 
(Lyzenga, 1978), referred to as photobathymetry in the table, is good for relatively clear 
water during cloudless, daytime passes and has one of the best horizontal resolutions, but 


17 






it is not a good option for turbid waters, cloudy regions, or nighttime passes. It also 
requires in situ knowledge for calibration purposes, negating its usefulness in inaccessible 
areas (Abileah, 2006). 

The wave refraction SAR method (Wackerman & Clemente-Colon, 2004) is a 
good all-weather, day or night option, and it does not rely on water clarity; but it has the 
least horizontal resolution and is not the most accurate method. The WKB method in this 
study, works fine in turbid waters, but is not a good option when the scene is obscured by 
clouds or for nighttime passes (when optical imagery is used), or when ocean waves are 
not present (harbors, calm seas, etc.). It has decent depth accuracy and horizontal 
resolution, and does not need in situ data for calibrating, making it well suited for denied 
areas. In order to optimize the wave contrast in the images, the camera should point in 
the direction of the sun while avoiding the 20-30° sun glint cone (Abileah, 2006). 

This study assesses a WKB computer algorithm that automates the nearshore 
bathymetry extraction on imagery of the southeastern portion of the Hawaiian island of 
Oahu, obtained by the WorldView-2 satellite. 

The S spectral bands of WorldView-2 

Pan 
Coastal 
Blue 
Green 
Yeltwr 
Red 
Red Edge 
NIRi 
NIRz 

3 SO 4 SO 550 650 750 £50 S50i tOSO 

Wave length (rim) 

Figure 8. WorldView-2 sensor bands showing their relative positions and 
overlap in the electromagnetic spectrum (From DigitalGlobe, 2011b). 
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MATERIALS 


1. WorldView-2 Satellite 

The WorldView-2 satellite was launched into a 770 km high, sun synchronous 
orbit in October 2009 from Vandenberg Air Force Base. The third in DigitalGlobe’s 
commercial imaging satellite constellation, it has a state-of-the-art multispectral optical 
imager and provides approximately 2 m multispectral and 0.5 m panchromatic resolution 
imagery. The panchromatic and eight multispectral bands and their relative positions and 
overlap in the electromagnetic spectrum are illustrated in Figure 8. The sensor has a 
revisit time of 1-4 days, depending on desired viewing angle and resolution. It collects 
imagery in a swath 16.4 km wide at nadir and can slew 200 km in 10 seconds. Some 
additional key specifications are listed in Table 2 (DigitalGlobe, 2011a). 


Table 2. Key WorldView-2 specifications (After DigitalGlobe, 2011a) 


Sensor Bands 

Panchromatic: 450 - 800 nm 

8 Multispectral: 

Coastal: 400 - 450 nm Red: 630 -690 nm 

Blue: 450 - 510 nm Red Edge: 705 - 745 nm 

Green: 510 - 580 nm Near-IRl: 770 - 895 nm 

Yellow: 585 - 625 nm Near-IR2 : 860 - 1040 nm 

Sensor Resolution 

Panchromatic: 0.46 m GSD at nadir, 0.52 m GSD at 20° off-nadir 
Multispectral: 1.85 m GSD at nadir, 2.07 m GSD at 20° off-nadir 

Dynamic Range 

11-bits per pixel 

Swath Width 

16.4 km at nadir 

Attitude Determination 
and Control 

3-axis stabilized 

Actuators: Control Moment Gyros (CMGs) 

Sensors: Star trackers, solid state IRU, GPS 

Pointing Accuracy 
and Knowledge 

Accuracy: <500 m at image start and stop 

Knowledge: Supports geolocation accuracy below 

Retargeting Agility 

Acceleration: 1.43 deg/s/s 

Rate: 3.86 deg/s 

Time to Slew 200 km: 10 sec 

Onboard Storage 

2199 Gb solid state with EDAC 

Communications 

Image and Ancillary Data: 800 Mbps X-band 

Housekeeping: 4, 16 or 32 kbps real-time, 524 kbps stored, X-band 
Command: 2 or 64 kbps S-band 

Max Viewing Angle / 
Accessible Ground 

Swath 

Nominally -*-/-45° off-nadir = 1651 km wide swath 

Higher angles selectively available 

Max Contiguous 

Area Collected in 
a Single Pass 
(30® off-nadir angle) 

Mono: 138 x 112 km (8 strips) 

Stereo: 63 x 112 km (4 pairs) 

Revisit Frequency 
(at 40*N Latitude) 

1.1 days at 1 m GSD or less 

3.7 days at 20° off-nadir or less (0.52 m GSD) 

Geolocation Accuracy 
(CE.90) 

Specification of 5.0 m CE90 at less than 30° off-nadir, with 
predicted performance in the range of 4.6 to 10.7 m (15 
to 35 ft) CE90, excluding terrain and off-nadir effects 

With registration to GCPs in image: 2.0 m (6.6 ft) 


19 


2. SHOALS Bathymetry 

In order to properly assess the results obtained with the WKB algorithm, a 
comparison to ground truth depths was necessary. LIDAR bathymetry, obtained for the 
purpose of coral reef mapping by JALBTCX using the CHARTS system, is provided as 
ground truth. The survey was conducted in 2000 using an Optech, Inc., SHOALS-3000 
LIDAR instrument integrated with an Itres CASI-1500 hyperspectral imager. SHOALS 
capabilities meet U.S. Army Corps of Engineers Hydrographic Survey accuracy 
requirements for Class 1 surveys and the IHO nautical charting standards for Order 1. 
The positional accuracy was +/- 3 m in the horizontal and +/- 15 cm in the vertical with a 
resolution of 0.00000001 degrees in latitude and longitude and 0.1 m in depth 
(JALBTCX, 2011). 

For comparison, it is preferable to have recent data to ensure that the bathymetry 
has not evolved. The SHOALS data are the most recent available for the area of interest 
and was collected over 10 years before the satellite pass that produced the images 
analyzed in this study. This may be a source of error, depending on bathymetric 
evolution during this interval. Another potential source of comparison error is UTM 
coordinate in-accuracies, causing slight differences in the pixel locations between the two 
data sets. Another consideration is the bathymetry derived from the satellite data set 
includes tide effects, whereas the SHOALS data are referenced to mean lower low water 
(MLLW). The satellite data must be corrected to account for tide height. According to 
National Oceanic and Atmospheric Administration tide tables, the tide varied between 
0.06 m at low tide and 0.50 m at high tide on March 31, 2011, and was approximately 
0.37 m above MLLW at Waimanalo Beach at the time of the collection. 

3. WKB Algorithm 

The WKB algorithm is a patent-pending method for generating maps of nearshore 
depth and surface currents from a variety of imaging inputs from various platforms 
(Abileah, 2011). As of December 11, 2011, it is composed of almost 100 subroutines 
written in MathWorks’ MATLAB software and is a work in progress with several 
updates planned over the coming months. 


20 




Figure 9. STK snap shot of WorldView-2 collection pass. 


C. METHOD 

The imagery was collected at about 2200Z on March 31, 2011, by the 
WorldView2 satellite. A Satellite Tool Kit (STK) snap shot shows the geometry part 
way through the collection (Figure 9). The satellite is viewing a swath of the southeast 
tip of Oahu and the adjacent ocean from the west looking toward the east in the direction 
of the sun, which is the ideal geometry for high signal-to-noise ratio (SNR) as long as the 
sun glint cone is avoided. The fact that the waves are propagating east to west toward the 
sensor also improves SNR. The scene is mostly cloud-free over the water and has 
adequate ocean waves present, both necessary conditions for WKB. The Waimanalo 
Beach area, covered by zone 8 in the SFIOALS bathymetry (Figure 10), was chosen as 
the test area. 

The imagery consists of two sets of three pairs of panchromatic and multispectral 
images taken approximately 10 seconds apart, and was obtained from DigitalGlobe via 
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the National Geospatial-Intelligence Agency. The images were provided in Basic IB, 
ortho-ready format and geographic latitude-longitude coordinates. An example image 
from the multispectral set is shown in Figure 11. The images were orthorectified and 
converted to UTM coordinates for the WKB algorithm. A summary of the multispectral 
images is given in Table 3. To take advantage of their ability to increase wave contrast 
by reducing subsurface noise and atmospheric path radiance, the red and near-infrared 
bands (5 through 7) are used for the bathymetry extraction. The panchromatic images, 
which are not listed in the table, have the same metadata; however, the filenames contain 
a “P” rather than an “M.” 



Figure 10. SHOALS 2000 survey of Oahu (From JALBTCX, 2011). 

The images are pre-processed prior to entering the WKB algorithm by cropping to 
a user-defined size, loading the cropped images into a data cube, verifying the time 
interval between successive images in the set, registering the images to each other, 
conducting a deep water spectral analysis, and determining which pixels are water and 
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which are not. Then the program conducts a WKB extraction on the data bundle(s) using 
the specified tile size(s) (this determines the horizontal resolution of the resulting depth 
and current fields). 


Table 3. Summary of image files. 


Image Name 

First Line Time 

Mean Sat Az 

Mean Off-nadir 

View Angle 

90-mean Sat El 

meanSatEl 

Set 1 






11MAR31215059-M1BS-052517305060 

21:50:59 

311.80 

38.5 

44.2 

45.8 

11MAR31215109-M1BS-052517305010 

21:51:09 

306.10 

37.4 

42.8 

47.2 

11MAR31215119-M1BS-052517305030 

21:51:19 

300.00 

36.5 

41.8 

48.2 

Set 2 






11MAR31215155-M1BS-052517305050 

21:51:55 

276.30 

35.9 

41.0 

49.0 

11MAR31215204-M1BS-052517305040 

21:52:05 

270.10 

36.4 

41.7 

48.3 

11MAR31215214-M1BS-052517305020 

21:52:14 

264.20 

37.2 

42.7 

47.3 



Figure 11. Example WorldView-2 multispectral image. 
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1 . 


Pre-Processing 


The pre-processing ensures the images are ready for the WKB algorithm. The 
current version is described, but it need not be done exactly as it is today. The 
multispectral and panchromatic image sets are processed separately. First, the image files 
in each set are inventoried. At least two are required and can be multispectral or 
panchromatic (but they have to be the same type with the current version). 

A graphical user interface (GUI) allows several user inputs to be entered during 
the process. The area for processing is defined by cropping a Google Earth picture 
(Figure 12). The cropped images are then loaded into analysis data cubes (one for each 
set) at a down-sampled resolution that increases efficiency but still resolves the smallest 
waves so as to not lose any wave energy (3 m for the multispectral and 1.5 m for the 
panchromatic data). During the loading process, the images are stitched together to 
create full scenes if they were partitioned into smaller strips by the vendor (such was the 
case with this imagery). 



6.25 


6.35 6.4 

UTM Easting 


6.45 


xIO 


Figure 12. Example cropped area (in green). 
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To get the timing precisely correct for the wave propagation, the time each scan 
line was imaged is detennined, based on the first line time and the scan rate in the image 
metadata. Without accurate timing, wave speed will not be correct and the resulting 
calculations for both current and depth will be off. 

Image registration is performed using fixed land features to properly align the 
image frames with each other. The area chosen should be close to sea level and have 
buildings or other fixed structures in it (Figure 13). Registration calculates the mis¬ 
registration and shifts the images to ensure the pixels at the same geographic location on 
each image are overlaid on top of the each other, to within sub-pixel accuracy. The 
panchromatic images are used to compute the registration shift because of their higher 
resolution, but the multispectral could also be used if panchromatic imagery is not 
available. The multispectral images are shifted by the same amount as the panchromatic 
ones. It should be noted that, if the registration is off slightly, it will simply lead to a 
current bias that does not impact the bathymetry detennination. 

Use mouse to outline registration chip... 

4000 

3500 
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2500 

2000 

1500 

1000 

500 

Figure 13. Registration GUI 
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A deep water box is then analyzed to determine the time interval between images 
and to compare the spectra of the multispectral and panchromatic sets (Figure 14). The 
interval is calculated using the fact that, per the dispersion relation as applied to deep 
water, wave speed is constant for a given frequency. Therefore, it is important to pick a 
spot where depth is at least half of a wavelength of the primary waves to get an accurate 
time interval (usually about 50 m). This calculation would be unnecessary if the timing 
metadata were more precise. It is important to determine the timing very precisely for the 
wave speed calculation. The interval timing is combined with the time offset determined 
earlier to obtain very accurate time spacing for pixels in different images. Other analyses 
are conducted with regard to radiance that are not important for this work, but may be of 
interest to researchers using radiometric techniques. 


Select deep water location MSI deep water chip Pan deep water chip 

1900 

1850 
1800 
1750 
1700 
1650 
1600 


1700 1800 1900 2000 3400 3600 3800 4000 

Pixels Pixels 

MSI k-spectrum Pan k-spectrum 





Figure 14. Deep water analysis 


The final pre-processing function creates a mask that prevents pixels that are not 
water from being included in the spectral analysis (Figure 15). The various bands of the 
multispectral images are compared to determine if a pixel most likely contains land, 
water, clouds, a boat or some other interfering phenomenon. With the “bad” pixels 
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removed, and all the other adjustments made during pre-processing, the images are ready 
to go into the WKB algorithm for extraction of the bathymetry. 

2. WKB Extraction 

WKB is rooted in the fact that surface gravity waves decrease in speed as they 
move to shallower depths, according to the linear dispersion relation given by Equation 1. 
By capturing images of ocean waves at short, precisely known time intervals, their 
horizontal displacement can be observed and celerity calculated. Depth is inferred from 
the dispersion relation. WKB becomes less useful in the cases where depths approach the 
height of the waves, such as in the surf zone (due to non-linear effects that are not 
incorporated in the linear dispersion relation); and where depths are greater than half the 
wavelength of the longest detectable waves (the region where wave speed is independent 
of depth). 



Figure 15. True color image (left) with corresponding land mask (right) 

The key to extracting bathymetry from the information contained in the images is 
the Fourier transform. Most commonly, 3-D Fourier analysis is applied to transform 


27 



image data from two-dimensional space-time into wavenumber-frequency spectra. This 
technique requires long time series lengths, or dwell times -100 seconds (Dugan et ah, 
2001a). 

Two-dimensional Fourier analysis is more appropriate for satellite imagery, 
which, due to the rapidly changing geometry, has effective dwell times equal to the time 
between images. The 2-D analysis allows for these much shorter ~10 second dwell times 
without adversely increasing the error. This is because the depth accuracy does not 
depend on dwell time so the frequency resolution is not relevant (Abileah & Trizna, 

2010) . In this case, the 2-D data are transformed into 2-D wavenumber spectra. In both 
methods, the image is broken into 2-D tiles, the dimensions of which are driven by a 
compromise between maximizing spatial resolution and minimizing depth error (Abileah, 

2011) . Tiles that are too small provide better resolution at the expense of depth accuracy. 
Large tiles produce better accuracy, but resolution suffers since one depth is calculated 
per tile. 

The 2-D algorithm, a modification of which is applied in this study, applies a 
propagation kernel to the 2-D Fourier transforms of the N images {^(x^S^x),..., S N (x)} 

,F n (k) - 3 2 [,SJ . The propagation kernel moves the waves forward or backward in time 
and is defined as, 


0 ± (k x ,k y \d,u) = e 


±ico 0 (k x ,k y \d,u)T 


( 11 ) 


where d is depth, u is the surface current, r is the time interval between compared 
images, and k x and k y are the x and y components of the wavenumber. The sign in the 
exponent detennines whether the waves propagate forward or backward. The term 
co 0 (k x ,k v | d,u) in the exponent comes from the linear gravity wave dispersion relation in 
this form: 


co 0 (k \d,u) = | k | tanh(| k\d)-k ■u 


( 12 ) 
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which also accounts for modification of the waves by ocean currents. Using the 
propagation kernel, the Fourier transfonn at time n is related to the next at time n+m by 

F n+m = 0"'F (Abileah & Trizna, 2010). (13) 

Wave energy in the transform space tends to be distributed along a dispersion 
surface and the depth and current can be found by finding the best fit of the dispersion 
formula to this measured power spectrum. One way to find the best fit is by minimizing 
the difference between successive Fourier transfonns by tuning the propagation kernel to 
the correct depth and current. Mathematically, this is done by minimizing the objective 
function 

N—m 

argmin, , J A = £ W(k)^ I F„„ | 2 

, (14) 

where W ( k ) is a weighting filter that increases the accuracy of the solution by increasing 
the SNR. It does this by masking much of the background noise by keeping those 
wavenumbers where the wave energy is concentrated and eliminating those where it is 
not. Since waves are located where the energy adds coherently in wavenumber space 
when the images are propagated to the same time, this determination is done by seeing 
for which wavenumbers this is true. The objective function also includes a provision for 
more than two image pairs by including a summation over frame intervals, in, from n = 1 
to N-m (Abileah, 2011). 

A patent-pending technique, which is a modification of the 2-D algorithm, 
overcomes the accuracy-resolution compromise discussed above by adding a 2-D inverse 
Fourier transform operation, 3 j 1 , to the objective function (Abileah, 2011). By 

transfonning the data from wavenumber space back into spatial coordinates, a depth is 
obtained for every image pixel rather than for every tile. In fact, it eliminates the need to 
tile the images, although this has not yet been incorporated into the code. To expedite 
implementing the new method, the existing code was modified rather than generating 
new code, which is computationally inefficient but yields the same results. A more 
efficient revision is forthcoming. 
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Besides better spatial resolution and less computation time, another benefit of the 
new method is it is less susceptible to interference from non-water pixels because, unlike 
in wavenumber space, in spatial coordinates the waves can be separated spatially from 
the interfering phenomena. 

As before, the 2-D Fourier transfonn is applied to the (in this case, tiled) images, 
converting them into 2-D wavenumber spectra. (New implementations can forgo the 
tiling step; however, tiling may still be desired to take advantage of the multiple 
processors in parallel computing settings.) Minimizing the revised objective function: 


ar § J a 0) = 3 2 ! \_ w ( k )( F i ( k ) - ° ' F 2( k )) 


(15) 


again by tuning ®, produces the best depth and current solution for each x (Abileah, 

2011 ). 

Although W eliminates much of the noise, some is still present after the filtering. 
To further improve the results, a sliding window A(x ) is applied. The window size can 
be adjusted as necessary to improve SNR; however, a larger window will reduce 
resolution, but not to the degree the tiles did in the old method. With the appropriate 
A(x ), spatial resolution with the new method can be 10 m (Abileah, 2011). The code 
implemented herein used 250 m tiles with 50% overlap, creating an effective spatial 
resolution of about 125 m. 

For instances where more than two images are available, SNR can be increased 
further by averaging across the images, n = 1 to N-m, not unlike in the old method. For 
data sets with long interval times, and as was done for this study with 10 seconds 
between images, consecutive images are used, meaning m= 1. Still another improvement 
through averaging is to sum over a set of interval spacing values, M, so multiple image 
separation options are available. The final, general form for the objective function that 
includes all these possibilities is (Abileah, 2011): 

argmin dJl J A (x)= £ £ £ [W{k){F, (k)-®-\F 2 (k))]| . (16) 

me A/ n =1 A(x) 
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Choosing an m other than one creates comparisons of images that are not 
consecutive. Since this increases the time interval between compared images, effectively 
decreasing the sampling rate, it should be done with the following consideration in mind. 
In wavenumber space, a phase shift is equivalent to wave displacement in spatial 
coordinates. Using Equation 11, the measured phase shift leads to the depth result 
obtained with Equation 16. In theory, a very small time interval is all that is necessary to 
determine the phase shift. In the real ocean, where noise complicates the wave signal, 
however, a longer time is required to accurately measure the phase shift amidst the 
background noise. On the other hand, too long of an interval allows the waves to become 
less coherent as they become modified by wind, refraction, bottom friction, etc. 
Therefore, a time interval of at least one second, but not more than several seconds, is a 
good compromise. Ideally, the 10-second interval time between the images in this study 
should be shorter to ensure less wave modification. With this in mind, an m larger than 
one would only be appropriate for data sets having image interval times (or sampling 
times) of one second or less between images. This can also be extended to using sets of 
different values of m (Abileah, 2011). 

There are a couple of ways to implement the objective function tuning, which 
involves searching for the depths and current speeds that minimize Equation 16. One 
method attempts this by trying every possible combination of depth and current speed. 
This option is computationally intensive since it has to loop through three variables, 
depth and both the x and y components of the current velocity, in computing the 
propagation kernels for use by Equation 16 (Abileah, 2011). 

Another way to do it that reduces the computation time significantly (and the way 
the current revision accomplishes it) is to separate the current search from the depth 
search. To do this, a low wavenumber filter is applied to isolate the higher wavenumbers 
(short wind waves) because they are effectively deep water waves and unaffected by 
depth. This removes the need to include depth in the first search and allows any 
deviation from the dispersion relation predicted wave propagation to be attributed to the 
current. Once the current is known, the images are shifted, similar to the registration 
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process, to remove any wave displacement caused by the current. Then the search for 
depth is done to complete the solution (Abileah, 2011). 

Of note, in execution of Equation 13, the current code does not propagate the 
waves in one image to the time of the other image being compared. Rather, it propagates 
the waves in both images, one forward and one backward, to the time half-way between 
the two (Abileah, 2011). This was done purely for cosmetic reasons. 

The output of the WKB program includes several figures, some of which are 
shown above, a diary that contains a chronological account of the program’s operations, 
image registration movies that show the shift being applied to the images, and image data 
cubes that contain the image data at various stages of the process. The main output figure 
is a quad chart that displays wave direction on the top left, a true color multispectral 
image of the area processed on the top right, the extracted bathymetry on the bottom left, 
and the extracted ocean currents on the bottom right (Figure 16). Another key figure, 
which will be discussed in Chapter IV, shows how the WKB derived bathymetry 
compares to the SHOALS bathymetry. 
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Figure 16. Image Set 1 case 3MSI WKB output showing wave direction (top 
left), a true color image (top right), extracted bathymetry (bottom left), and 
extracted ocean currents (bottom right) 
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IV. ANALYSIS RESULTS 


A. GENERAL APPROACH 

For the comparison analysis of WKB extracted bathymetry to SHOALS 
bathymetry, two data runs were conducted by ingesting two different image sets into the 
WKB algorithm. The first set consisted of the first three multispectral images in Table 3 
and their corresponding panchromatic images, while the second set consisted of the last 
three multispectral images and their corresponding panchromatic images. While it was 
not practical to make the user-selectable items identical between runs, the cropped area of 
interest, registration box, and deep water box were carefully chosen to minimize 
differences between the two sets. This way, any differences in the results would be 
mostly attributable to differences in the images. 

Four cases were produced in each of the two data runs, two using multispectral 
images and two using panchromatic images. The four cases, named for how many of 
each type of image were ingested, are: two multispectral images (2MS1), three 
multispectral images (3MS1), two panchromatic images (2Pan), and three panchromatic 
images (3Pan). For the two-image cases, the first two images in each set were the ones 
ingested. 

The WKB extracted depth fields were compared to the SHOALS bathymetry 
using side-by-side area maps. These figures show the SHOALS map on the left and the 
WKB map on the right and allow a quick visual comparison of the WKB field with 
ground truth field (Figure 17). 

East-west transects were chosen to investigate a few sample cross-sections for 
each case. For each transect, the WKB depth is plotted with the SHOALS depth as a 
measure of the similarity between them along the east-west slice. The depth error is 
calculated for each transect, as well. This graph aids in identifying where along the 
transect the greatest deltas are located. A third graph for each transect is the relative 
depth error, which shows depth error as a percentage of depth. This allows visualization 
of where the error has more or less significance (Figure 18). 
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Figure 17. Comparison maps of SHOALS (left) and 3MSI WKB for Image 

Set 1 (right) 


In addition to transects, each case was analyzed using all the valid depths in the 
entire subject area. Scatter plots show the distribution of the data and their correlation. 
Since the data field consists of millions of points, resulting in saturated plots, the data 
were thinned out by picking out every thousandth point to plot. The plots specify the 
slope of the linear regression line, its y-intercept, and R , the square of the correlation 
coefficient for the SHOALS-observed and WKB-predicted values. 

The regression line slope and R 2 are indicators of WKB performance. If WKB 
were perfect, its depths would be exactly correlated to the SHOALS depths, which would 
yield slope and R values of 1.0. The less perfectly WKB performs, the less correlated 
the WKB and SHOALS depths will be, and the further these values decrease from 1.0. 
The scatter plots have a second panel that filters out all SHOALS depths greater than 15 
m, for reasons discussed later (Figure 19). 
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Figure 18. Transect from Image Set 1 case 3MSI showing WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 19. Scatter plot from Image Set 1 case 3MSI showing thinned data for 
all depths (top) and for just depths less than 15 m (bottom) 
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In addition, a bar graph was produced that shows, for each case, the mean error 
for several depth partitions, or bins. The bins are one to 2 m, followed by every 2 m up 
to 20 m, with the last bin including all depths greater than 20 m. This aids in determining 
at which depths the algorithm worked better (Figure 20). 


WKB Depth Error vs Depth 



Depth (m) 

Figure 20. Bar graph from Image Set 1 case 3MSI showing the mean depth 

error for several depth bins 


B. IMAGE SET 1 

A quick glance at the comparison map for 3 MSI (Figure 17) shows that, overall, 
WKB did a good job of depicting the main bathymetric features, such as the elongated 
trough (yellow feature) near the center of the map and the shallower spots to the south of 
it (dark red). The abrupt, steeper slope beyond the more gently sloping beach is also well 
depicted where the red quickly turns yellow a couple of kilometers from the beach. 
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Image Set 1 results are summarized in a line graph that consolidates the mean 
depth errors in the bar graphs of all four cases on one chart (Figure 21), in a table listing 
the mean depth error and relative error ranges for each depth bin (Table 4), and in a table 
listing the R 2 values for each case (Table 5). The results highlight some of the strengths 
and weaknesses of the method. The complete series of graphs for Image Set 1 is found in 
Appendix A. 


WKB Depth Error vs Depth 



Figure 21. Mean depth errors for all four Image Set 1 cases 

Understandably, the 1-2 m bin for all cases is dominated by extremely large mean 
depth errors. The nonlinear processes in the surf zone render the linear dispersion 
relation inaccurate, and therefore, WKB less effective at predicting depths below about 
2 m. To improve WKB in this regime would require it to account for nonlinearities as the 
depth approaches the height of the waves. This could be done by applying a correction 
within the linear dispersion relation (Fledges, 1976), but would still be an approximation. 
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Other factors contributed to error at these shallow depths, including an ambiguity caused 
by a sampling time that is equivalent to the period of some of the waves, and insufficient 
short waves for current estimation. 


Table 4. Mean depth error and relative error ranges for each depth bin in Image Set 1 


Case 

<2m 


2-4m 



4-6m 



6-8m 


2MSI 

6.8 

341% 

3.8 

94% 

188% 

0.9 

14% 

21% 

1.0 

13% 17% 

3MSI 

8.4 

420% 

1.6 

40% 

80% 

0.6 

11% 

16% 

0.8 

10% 14% 

2Pan 

3.9 

196% 

1.7 

42% 

84% 

1.3 

22% 

33% 

1.9 

24% 32% 

3Pan 

5.8 

288% 

1.1 

28% 

55% 

1.2 

20% 

29% 

1.3 

16% 22% 












Case 

8-10m 



10-12m 



12-14m 




2MSI 

1.4 

14% 

17% 

2.3 

19% 

23% 

3.6 

26% 

30% 


3MSI 

1.5 

15% 

18% 

2.5 

21% 

25% 

2.1 

15% 

18% 


2Pan 

3.0 

30% 

37% 

2.4 

20% 

24% 

3.4 

24% 

29% 


3Pan 

2.9 

29% 

37% 

3.1 

26% 

31% 

1.9 

14% 

16% 



Case 

14-16m 



16-18m 



18-20m 



>20m 


2MSI 

4.7 

29% 

33% 

4.5 

25% 

28% 

5.6 

28% 

31% 

4.7 

23% 

3MSI 

3.8 

24% 

27% 

4.3 

24% 

27% 

5.6 

28% 

31% 

4.8 

24% 

2Pan 

4.6 

29% 

33% 

6.6 

37% 

41% 

6.6 

33% 

37% 

5.2 

26% 

3Pan 

3.5 

22% 

25% 

4.8 

27% 

30% 

6.1 

31% 

34% 

4.8 

24% 


2 

Table 5. R~ values for Image Set 1 


Case 

R 2 

R 2 < 15 m 

2MSI 

0.8681 

0.8909 

3MSI 

0.9005 

0.9078 

2Pan 

0.8603 

0.8453 

3Pan 

0.905 

0.911 


40 



Another limitation occurs at the other end of the depth field, where mean depth 
errors were relatively high. At depths greater than half the wavelength, waves are not 
affected significantly by the bottom, so they are dispersive and the dispersion relation is 
not helpful for determining depth. Therefore, WKB can only determine depth up to the 
point where the water depth allows gravity waves to become nondispersive. If fewer 
longer waves are present in the wave field, this will occur at shallower depths. If longer 
waves are abundant, then WKB will be effective at deeper depths. 

Therefore, WKB is limited on its ability to extract depth as determined by the 
wave spectrum detected by the sensor. In this instance, it appears that the spectrum lacks 
sufficient energy in waves longer than about 40 m, so that by a depth of about 16 m, the 
mean depth errors are consistently above 4 m with relative errors of 25% and greater. 
This assumption is supported by The Coastal Data fnfonnation Program (CD1P) of 
Scripps Institution of Oceanography, which indicates the sea on 31 March off of nearby 
Mokapu Point was relatively calm with very little long wave swells (Figures 22 and 23). 
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Day of Month (UTC) 


098 MOKAPU POINT, HI 

(BUOY) 

ENERGY SPECTRUM 
MARCH 2011 



rttp/ZcdipLcadedu 1 ’ 


Figure 22. Energy spectrum from Mokapu Point buoy for March 2011, 
showing larger swell events earlier in the month, but relatively little swell on 
the date of collection (31 March) (From CDIP, 2012) 
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Figure 23. Energy spectra from Mokapu Point buoy for the time of collection 
(left) and a swell event a few days prior (right), showing an order of 
magnitude difference in energy at the low frequencies (From CDIP, 2012) 

In between these limits, perfonnance was much better. 3MSI had the best overall 
performance, with a mean depth error of 1.8 m from 2-16 m, a regression line slope of 
1.01 and R 2 of 91%. The next best was 3Pan with 2.1 m mean depth error, slope of 1.06 
and R " of 91%. The two-image cases did the worst with mean depth errors of 2.5 m and 
2.6 m for 2MSI and 2Pan, respectively. It was expected that the three-image cases would 
perform better than the two-image ones because of the averaging technique that increases 
SNR with multiple image pairs. However, the multispectral outperfonning the 
panchromatic was unexpected and is addressed shortly. 

In the shallower depths outside the surf zone, the multispectral cases performed 
quite well. From about four to 10 m, mean depth errors were 1.5 m, or less, and relative 
errors were in the 10-20% range. By 10 m, the best panchromatic case (2Pan in this 
depth bin) did as well as the multispectral cases, all having mean depth errors of about 
2.4 m and relative errors in the 20-25% range. A bit of an anomaly, this was the only 
instance where 2Pan was better than 3Pan. Looking at the comparison maps in Appendix 
A, it may be due to 3Pan doing a slightly worse job in the area west of Manana Island in 
the southeast corner, an area which is primarily 10-12 m in depth. 
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The expectation was that the panchromatic cases would outperfonn the 
multispectral cases because of their higher resolution images being able to better resolve 
the shorter waves leading to a more accurate current determination. This, in turn, would 
produce more accurate depths. A possible reason for the multispectral cases’ better 
performance has to do with water clarity. The ocean bottom can be seen quite clearly in 
the shallower regions in Figure 11. This certainly added to the non-gravity wave noise in 
the panchromatic cases with a corresponding reduction in SNR. The multispectral cases 
were partially immune to this problem by screening out most of the subsurface noise by 
using only bands 5 through 7. What the multispectral cases lacked in resolution, they 
made up for in band selectivity and vice versa for the panchromatic cases. In the deeper 
water, where bottom radiance was less significant, the panchromatic cases did better. 

In line with this reasoning, 3Pan perfonned the best in the 12-16 m region, edging 
out 3MSI by a couple tenths of a meter in depth error. Mean depth errors were about 2 m 
in the 12-14 m bin and 3.5 m in the 14-16 m bin for these two cases. The worse¬ 
performing two-image cases had well over 3 m errors in the 12-14 m bin and well over 4 
m in the 14-16 m bin. 

C. IMAGE SET 2 

Image Set 2 yielded similar results to set one, with some notable differences. 
Again, WKB did a good job of depicting the main bathymetric features, as seen in the 
comparison maps in Appendix B, where a complete series of Image Set 2 graphs are 
found. As was done for Image Set 1, Image Set 2 results are summarized in a line graph 
of all four cases that consolidates the mean depth errors on one chart (Figure 24), in a 
table listing the mean depth error and relative error ranges for each depth bin (Table 6), 
and in a table listing the R~ values for each case (Table 7). 

The same limitation in the surf zone is apparent in Image Set 2’s results, with the 
1-2 m bin for all cases dominated by extremely large mean depth errors. In the deeper 
region, although the same limitations applied to this set, better results were achieved with 
set two. For instance, a respectable mean depth error of 3.0 m and corresponding relative 
error of 17-19% was obtained by 3MSI in the 16-18 m bin. In comparison, all cases in 
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set one had errors well above 4.0 m and greater than 23% by 16 m depth. In fact, 3MSI 
performed better in the 18-20 m bin in set two than it did in the 16-18 m bin in set one. 
It was not until the final greater-than-20 m bin that 3MSI’s performance degraded above 
the poorer set one values. The performance may have been aided by geometry more 
favorable for enhancing long wave contrast. If this was the case, the wave spectrum has a 
higher SNR in the long wave portion, allowing for deeper bathymetry extraction. 
However, this is purely speculative. 


WKB Depth Error vs Depth 



Figure 24. Mean depth errors for all four Image Set 2 cases 

As in set one, the three-image cases outperformed their two-image counterparts 

almost exclusively. Also like set one, 3MSI had the best overall performance in between 

the surf zone and the deeper limit, with a mean depth error of 1.7 m from 2-16 m, and an 

only slightly larger 1.9 m when expanded to include 16-18 m. Regression line slope and 

R~ for 3MSI were 1.13 and 88%, respectively. The next best for set two was 2MSI, 

rather than 3Pan as was the case with Image Set 1. It had a mean depth error of 2.0 m for 
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the range of 2-16 m and similar regression line slope and R . The panchromatic cases did 
the worst this time with 2.6 m and 3.4 m for 3Pan and 2Pan, respectively. These figures 
are quite a bit worse than the 2.1m and 2.6 m achieved with set one. 


Table 6. Mean depth error and relative error ranges for each depth bin in Image Set 2 


Case 

<2m 


2-4m 



4-6m 



6-8m 



2MSI 

7.4 

369% 

1.6 

40% 

80% 

0.8 

13% 

20% 

1.0 

13% 

17% 

3MSI 

5.1 

256% 

1.6 

41% 

82% 

0.7 

11% 

16% 

1.4 

18% 

23% 

2Pan 

18.1 

903% 

4.4 

111% 

222% 

2.8 

47% 

71% 

3.4 

43% 

57% 

3Pan 

14.5 

726% 

3.1 

77% 

155% 

3.2 

54% 

80% 

2.7 

34% 

45% 













Case 

8-10m 



10-12m 


12-14m 




2MSI 

2.6 

26% 

33% 

2.3 

19% 

23% 

2.0 

14% 

17% 



3MSI 

2.2 

22% 

27% 

1.8 

15% 

18% 

1.8 

13% 

15% 



2Pan 

3.8 

38% 

48% 

2.9 

24% 

29% 

3.5 

25% 

29% 



3Pan 

3.6 

36% 

45% 

1.7 

14% 

17% 

1.5 

11% 

13% 


























Case 

14-16m 



16-18m 


18-20m 


>20m 


2MSI 

3.4 

21% 

24% 

3.8 

21% 

24% 

4.6 

23% 

25% 

5.4 

27% 

3MSI 

2.4 

15% 

17% 

3.0 

17% 

19% 

3.9 

20% 

22% 

5.0 

25% 

2Pan 

3.2 

20% 

23% 

5.3 

29% 

33% 

5.4 

27% 

30% 

5.0 

25% 

3Pan 

2.6 

16% 

19% 

4.7 

26% 

30% 

4.9 

25% 

27% 

4.3 

22% 


Table 7. R 2 values for Image Set 2 


Case 

R 2 

R 2 < 15 m 

2MSI 

0.8825 

0.8858 

3MSI 

0.9092 

0.8825 

2Pan 

0.8096 

0.8595 

3Pan 

0.8627 

0.8683 
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After examining the multispectral true color images on the top right of the WKB 
output figures in Appendices A and B, it is apparent that the ocean bottom appears 
brighter in the Image Set 2 images. This is most likely due to the difference in collection 
angle as the spacecraft moved over the area, with the angle becoming more favorable for 
bottom reflection as the pass evolved. This would have exacerbated the bottom generated 
noise issue for the panchromatic cases, leading to the further degradation in performance. 

The individual depth bin values tell a similar story, with the multispectral cases 
doing considerably better than the panchromatic cases through a depth of 10 m, where 
mean depth errors were nominally 1-2 m better and relative errors were 10-40% better in 
the multispectral cases. Above 10 m, the 3Pan case performed better than the 3MSI case, 
with mean depth errors of 1.5-1.7 m, until 14 m depth where 3MSI took honors again. 
Strangely, the multispectral cases did not degrade as rapidly as the panchromatic cases 
above 16 m, unlike in set one where they degraded roughly equally. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A WKB algorithm was investigated to determine its effectiveness in extracting 
bathymetry from WorldView-2 satellite optical imagery. This initial analysis indicates 
that the technique described herein is promising. The analysis involved comparing WKB 
extracted depths to SHOALS LIDAR survey obtained depths at Waimanalo Beach in 
Hawaii. 

The algorithm performed well in the midrange depths of 2-14 m, its best cases 
managing as little as 0.6 m and no worse than 2.3 m of mean depth error in any depth bin. 
Excluding the 2-4 m bin, relative errors for the best cases ranged from a best of 10% to a 
worst of 27%, with most being in the 10-17% range. The best regression line slope 
achieved was 1.01 with an associated R~ of 91%. 

The best analysis results came from the use of three multispectral images, the 
worst results came from the use of two panchromatic images. The panchromatic cases 
would certainly have performed better in a more turbid water or dark bottom area, where 
the bottom would not have contributed to a poor wave SNR. The three-image cases did 
better than the two-image cases in all but a couple of instances. 

There may have been some error introduced into the comparisons through a 
conversion of the original imagery from latitude-longitude geographic coordinates to 
UTM coordinates, which may have differed from the SHOALS UTM coordinates by 
several meters. Understanding that the bathymetry hardly changes in that amount of 
horizontal distance, and that the WKB algorithm produced depths at an effective spatial 
resolution of about 125 m, this error should be insignificant. Another, perhaps more 
significant, source of error is the elapsed time between the SHOALS survey and the 
imagery collection. During this almost 11-year period, the bathymetry may have 
changed, perhaps significantly in some portions. 

In spite of these potential sources of error and given the fairly complex 
environment in which it was tested, the performance of the algorithm was notable. 

49 



Unlike many previous WKB studies, which were conducted in areas with simpler, nearly- 
constant sloping bottoms, often assumed cross-shore currents, and avoided wave 
refraction, Waimanalo Beach has fairly complex bottom topography, intricate surface 
currents, and noticeable wave refraction. 

WKB is limited in its ability to do depth extraction in the surf zone (less than 
about 2 m) and beyond some deep limit determined by the surface wave spectrum. In 
this case, the lack of long wave swells limited WKB’s ability to determine depth 
satisfactorily to approximately 15 m. A scene with more long wave swells in it probably 
would have produced more accurate depths out to 20 m or more. 

B. RECOMMENDATIONS FOR FUTURE WORK 

The WKB algorithms and codes are still evolving. Continuing work will 
determine whether these changes improve the algorithm’s results. Among the newly 
implemented updates is the fusion of panchromatic and multispectral images into the 
same depth extraction routine, whereas they were calculated separately before. This 
allows the algorithm to take advantage of the strengths of each image set, using the 
higher resolution panchromatic images for the current determination and the bands of the 
multispectral images that create the best SNR to determine the depth. 

The code based on the tiling approach, used here, was recently replaced with code 
that makes the algorithm more computationally efficient with the non-tiling method. One 
way to really capitalize on this is to introduce parallel processing into the process. This 
would significantly speed up the computing time for larger areas and more image pairs. 

Future improvements include: a wave contrast metric to help identify when 
contrast has played a role in the results; better logic to improve performance with 
multiple image pairs; a way to flag extremes and replace them with appropriate values to 
eliminate unreasonable depth conclusions in the surf zone or at the deep limits; and 
accuracy and spatial resolution improvements. 

In addition to testing algorithm improvements, future work could investigate 
different combinations of multispectral bands to see which produces the best results by 
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better isolating the waves from the background. Perhaps principle component transforms 
of the images can be applied to increase wave visibility and SNR. 

Improvements in data collection could also be explored. For example, reducing 
the time interval between subsequent image pairs would reduce wave propagation and the 
change in view angle from one image to the next. This will help with accuracy by 
reducing the opportunity for the waves to be modified by sources other than depth 
(making the dispersion relation a more exact predictor of the wave motion) and by 
allowing the multiple-image-pair averaging technique to have more success. It would 
also alleviate the ambiguity issue associated with sampling at the same time interval as 
the period of a bulk of the waves. 

There are other ways to realize the last two gains if reducing the image collection 
time interval becomes too difficult. By incorporating a view-angle transfer function to 
account for the change in angle, the averaging technique could still produce improvement 
when averaging the image pairs. The ambiguity problem can probably be overcome with 
a better algorithm to more effectively deal with it. 

Finally, better choices in the imagery product could improve future efforts. 
Receiving the imagery in the UTM format that WKB requires would eliminate the need 
to convert it, avoiding any error introduced in the process. Obtaining accurate timing 
details from the vendor to establish more precisely when each image line was scanned 
will remove the need to estimate during the pre-processing. While the estimation is good, 
knowing the timing exactly will help improve WKB accuracy. 
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Figure 25. 2MSI WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 26. Comparison maps of SFIOALS (left) and 2MSI WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 27. 2MSI scatter plot showing thinned data for all depths (top) and for 

just depths less than 15 m (bottom) 
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WKB Depth Error vs Depth 
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Figure 28. 2MSI bar graph showing the mean depth error for several depth 

bins 


55 
























Relative Depth Error (%) Depth Error (m) Depth (m) Relative Depth Error (%) Depth Error (m) Depth (m) 


WKB-SHOALS Depth Comparison at 2360 UTM Northing (km) 





Figure 29. Transect at specified location showing 2MSI WKB (see below for 

more explanation of the panels) 
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Figure 30. Transect at specified location showing 2MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 31. Transect at specified location showing 2MSI WKB (see below for 

more explanation of the panels) 
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Figure 32. Transect at specified location showing 2MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 33. 3MSI WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 34. Comparison maps of SHOALS (left) and 3MSI WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



3MSI scatter plot showing thinned data for all depths (top) and for 
just depths less than 15 m (bottom) 
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Figure 35. 









WKB Depth Error vs Depth 
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Figure 36. 3MSI bar graph showing the mean depth error for several depth 

bins 


60 























Relative Depth Error (%) Depth Error (m) Depth (m) Relative Depth Error (%) Depth Error ( m ) Depth ( m ) 


WKB-SHOALS Depth Comparison at 2360 UTM Northing (km) 





Figure 37. Transect at specified location showing 3MSI WKB (see below for 

more explanation of the panels) 
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Figure 38. Transect at specified location showing 3MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 39. Transect at specified location showing 3MSI WKB (see below for 

more explanation of the panels) 


WKB-SHOALS Depth Comparison at 2363 UTM Northing (km) 
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Figure 40. Transect at specified location showing 3MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 41. 2Pan WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 42. Comparison maps of SHOALS (left) and 2Pan WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 43. 2Pan scatter plot showing thinned data for all depths (top) and for 

just depths less than 15 m (bottom) 
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WKB Depth Error vs Depth 
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Figure 44. 2Pan bar graph showing the mean depth error for several depth 

bins 
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Figure 45. Transect at specified location showing 2Pan WKB (see below for 

more explanation of the panels) 
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Figure 46. Transect at specified location showing 2Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 47. Transect at specified location showing 2Pan WKB (see below for 

more explanation of the panels) 
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Figure 48. Transect at specified location showing 2Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 49. 3Pan WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 50. Comparison maps of SHOALS (left) and 3Pan WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 51. 3Pan scatter plot showing thinned data for all depths (top) and for 

just depths less than 15 m (bottom) 
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WKB Depth Error vs Depth 
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Figure 52. 3Pan bar graph showing the mean depth error for several depth 

bins 
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Figure 53. Transect at specified location showing 3Pan WKB (see below for 

more explanation of the panels) 
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Figure 54. Transect at specified location showing 3Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 55. Transect at specified location showing 3Pan WKB (see below for 

more explanation of the panels) 


WKB-SHOALS Depth Comparison at 2363 UTM Northing (km) 
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Figure 56. Transect at specified location showing 3Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 57. 2MSI WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 58. Comparison maps of SFIOALS (left) and 2MSI WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



2MSI scatter plot showing thinned data for all depths (top) and for 
just depths less than 15 m (bottom) 
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Figure 59. 










WKB Depth Error vs Depth 
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Figure 60. 2MSI bar graph showing the mean depth error for several depth 

bins 
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Figure 61. Transect at specified location showing 2MSI WKB (see below for 

more explanation of the panels) 
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Figure 62. Transect at specified location showing 2MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 63. Transect at specified location showing 2MSI WKB (see below for 

more explanation of the panels) 
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Figure 64. Transect at specified location showing 2MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 65. 3MSI WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 66. Comparison maps of SFIOALS (left) and 3MSI WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



3MSI scatter plot showing thinned data for all depths (top) and for 
just depths less than 15 m (bottom) 
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Figure 67. 










WKB Depth Error vs Depth 
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Figure 68. 3MSI bar graph showing the mean depth error for several depth 

bins 
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Figure 69. Transect at specified location showing 3MSI WKB (see below for 

more explanation of the panels) 
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Figure 70. Transect at specified location showing 3MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 71. Transect at specified location showing 3MSI WKB (see below for 

more explanation of the panels) 
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Figure 72. Transect at specified location showing 3MSI WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 73. 2Pan WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 74. Comparison maps of SHOALS (left) and 2Pan WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 75. 2Pan scatter plot showing thinned data for all depths (top) and for 

just depths less than 15 m (bottom) 
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Figure 76. 2Pan bar graph showing the mean depth error for several depth 

bins 
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Figure 77. Transect at specified location showing 2Pan WKB (see below for 

more explanation of the panels) 
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Figure 78. Transect at specified location showing 2Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 79. Transect at specified location showing 2Pan WKB (see below for 

more explanation of the panels) 
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Figure 80. Transect at specified location showing 2Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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Figure 81. 3Pan WKB output showing wave direction (top left), a true color 
image (top right), extracted bathymetry (bottom left), and extracted ocean 

currents (bottom right) 
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Figure 82. Comparison maps of SHOALS (left) and 3Pan WKB (right) 


WKB-SHOALS Comparison 



WKB-SHOALS Comparison with a cutoff of-15 m 



Figure 83. 3Pan scatter plot showing thinned data for all depths (top) and for 

just depths less than 15 m (bottom) 
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Mean Depth Error (m) 


8 


WKB Depth Error vs Depth 



Figure 84. 3Pan bar graph showing the mean depth error for several depth 

bins 
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WKB-SHOALS Depth Comparison at 2360 UTM Northing (km) 





UTM Easting (km) 


Figure 85. Transect at specified location showing 3Pan WKB (see below for 

more explanation of the panels) 


WKB-SHOALS Depth Comparison at 2361 UTM Northing (km) 





UTM Easting (km) 


Figure 86. Transect at specified location showing 3Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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WKB-SHOALS Depth Comparison at 2362 UTM Northing (km) 





Figure 87. Transect at specified location showing 3Pan WKB (see below for 

more explanation of the panels) 


WKB-SHOALS Depth Comparison at 2363 UTM Northing (km) 





UTM Easting (km) 


Figure 88. Transect at specified location showing 3Pan WKB depth plotted 
with SHOALS depth (top); the difference between them, or depth error 
(mid); and depth error as a percentage of depth, or relative depth error (hot) 
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